# ---- Caption and Heading of File ---------------------------------------------
#                                                                              #
#        Replication script for analysis contained in:                         #
#        Measuring Fairness, by Barry C. Edwards                               #
#                                                                              #
#        Published in Alabama Law Review                                       #
#        Forthcoming: 2026                                                     #
#                                                                              #
################################################################################

# ---- Set-up machine, load datasets, helper functions, print settings      ----

  rm(list=ls())  # starting with clean environment

  # install the following packages: sate, RCPA3
  library(sate)  # make sure you're using version 2.3.0
  library(RCPA3)
  library(png)
  source("helperFunctions.R") 

  saveFiguresDesktop <- FALSE
  messageViewerFigures <- "Note: Figures look different in viewer compared to saved PNG files."
  verbose <- FALSE  # T/F switch to print more to console

    # set working directory to location of downloaded datasets on your machine
  setwd("C:/Users/Barry Edwards/Dropbox/jury research/Measuring Fairness Article") #  
  
  # Load dataset, already included in sate package (v. 2.3.0 and higher)
  all_juries_data <- observed.deliberations
  
# ---- Table 1: Verdict Preferences Estimated in Real Criminal Trials ----
  
  # reports P(g|actual) and P(g|hypo values in eight trials): 
  # U.S. v. King comes from Simon (1967) study
  king_verdict_preferences <- c("U.S. v. King", "Jury instructions", .760, .640)
  # Cal. v. Chapman from Winkelman et al (2014) study
  chapman_verdict_preferences <- c("Cal. v. Chapman", "Improper prosecutorial arguments", .600, .530)
  # Pa. v. Lee from Winkelman et al (2014) study
  lee_verdict_preferences <- c("U.S. v. King", "Improper prosecutorial arguments", .680, .660)
  # Pa. v. Berkowitz from Kahan (2010) study
  berkowitz_verdict_preferences <- c("U.S. v. Berkowitz", "Jury instructions", .650, .530)
  # Fla. v. Porter from my study
  surveyToAnalze <- "Porter" 
  source("analyzeOrginalSurveyData.R")
  porter_results <- survey_analysis_results
  porter_verdict_preferences <- c("Fla. v. Porter", "Ineffective assistance of counsel", 
                                  porter_results$compare_juror_table$actual_trial["P(g|actual)"], 
                                  porter_results$compare_juror_table$hypo_trial["P(g|hypo)"])
  # Fla. v. Washington from my study
  surveyToAnalze <- "Washington"
  source("analyzeOrginalSurveyData.R")
  washington_results <- survey_analysis_results
  washington_verdict_preferences <- c("Fla. v. Washington", "Ineffective assistance of counsel", 
                                      washington_results$compare_juror_table$actual_trial["P(g|actual)"], 
                                      washington_results$compare_juror_table$hypo_trial["P(g|hypo)"])
  # Ariz. v. Fulminante from my study
  surveyToAnalze <- "Fulminante"
  source("analyzeOrginalSurveyData.R")
  fulminante_results <- survey_analysis_results
  fulminante_verdict_preferences <- c("Ariz. v. Fulminante", "Coerced confession evidence", 
                                      fulminante_results$compare_juror_table$actual_trial["P(g|actual)"], 
                                      fulminante_results$compare_juror_table$hypo_trial["P(g|hypo)"])
  # Tex. v. Hopkins from my study
  surveyToAnalze <- "Hopkins"
  source("analyzeOrginalSurveyData.R")
  hopkins_results <- survey_analysis_results
  hopkins_verdict_preferences <- c("Tex. v. Hopkins", "Coerced confession evidence", 
                                   hopkins_results$compare_juror_table$actual_trial["P(g|actual)"], 
                                   hopkins_results$compare_juror_table$hypo_trial["P(g|hypo)"])

  verdictPreferencesTable <- data.frame(rbind(king_verdict_preferences,
                                              porter_verdict_preferences, washington_verdict_preferences,
                                              fulminante_verdict_preferences, hopkins_verdict_preferences,
                                              chapman_verdict_preferences, lee_verdict_preferences,
                                              berkowitz_verdict_preferences))
  colnames(verdictPreferencesTable) <- c("Study", "Error/Issue", "P(g|actual)", "P(g|hypo)")
  verdictPreferencesTable$`P(g|actual)` <- as.numeric(verdictPreferencesTable$`P(g|actual)`)
  verdictPreferencesTable$`P(g|hypo)`   <- as.numeric(verdictPreferencesTable$`P(g|hypo)`)
  print_table_neatly(verdictPreferencesTable)
  
  
# ---- Figure 1: Probabilities of Initial Guilty Vote Counts Given Jury Pool Preferences ----
figure1_info <- list(name="Binominal Theorem Illustration.png", width=6.5, height=2.25)
if(saveFiguresDesktop) { prettyPNG(figure1_info) } else { message(messageViewerFigures) }
  
    par(mar=c(3, 3, 2, 1), mfrow=c(1,2), family="serif", mgp=c(2, .65, 0), omi=c(0,0,0,0))
    binomialDistPlot(pg=.76, jury_n=12) # left panel
    binomialDistPlot(pg=.64, jury_n=12) # right panel
  
if(saveFiguresDesktop) { dev.off() }


# ---- Figure 2: Illustration of Deliberation between Factions ----
# This figure is not computed as it is a general illustration of concept.
figure2_info <- list(name="tug of war model.png", width=6, height=2)   
if(saveFiguresDesktop) { prettyPNG(figure2_info) } else { message(messageViewerFigures) }
    
  source("createTugOfWarIllustration.R")
  
if(saveFiguresDesktop) { dev.off() }
    
  
# ---- Figure 3: Jury Deliberation Process as a Markov Chain ----
figure3_info <- list(name="Markov Chain Illustration any size replicate.png", width=6, height=2)   
if(saveFiguresDesktop) { prettyPNG(figure3_info) } else { message(messageViewerFigures) }
  
  source("createGeneralMarkovChainDiagram.R")
  
if(saveFiguresDesktop) { dev.off() }
  


# ---- Table 2. Probability of Guilty Verdict Based on Jury’s Initial Poll ----
# the P(G|k) values are solved by a function that implements method described in article 
  print_table_neatly(data.frame(init_vote=0:12, 
             PG_n12=get_pG_by_k(jury_n = 12), 
             PG_n6=c(get_pG_by_k(jury_n = 6), rep(NA, 6))))


# ---- Table 3. Comparison of Empirical and Formal Models of Jury Deliberations ----

# Initial Model, Empirical Data
#  logregC(guilty_verdict ~ prop_jurors_g, data=all_juries_data, fit.stats=T, pre=T)
initial_empirical_model = glm(guilty_verdict ~ prop_jurors_g, family="binomial", data=subset(all_juries_data, !is.na(guilty_verdict)))

# Controlled Analysis, Empirical Data
z_empirical_model = glm(guilty_verdict ~ prop_jurors_g + death_penalty + six_person_jury, family="binomial", data=subset(all_juries_data, !is.na(guilty_verdict)))

# because formal model outcomes are Markov process, there is randomness. therefore
# the formal model is estimated many time and results (coef and fit stats) are averaged.
# storage objects created, filled by loops, and then averaged
  formal_estimates <- array(NA, dim = c(2,4,100))
  formal_fit_stats <- matrix(NA, nrow = 100, ncol = 4)
  # z_ variables are for the formal model with control variables
  z_formal_estimates <- array(NA, dim = c(4,4,100))
  z_formal_fit_stats <- matrix(NA, nrow = 100, ncol = 4)
  set.seed(12345)
  # the next block takes a while, depending on your processing speed
  for(j in 1:100) # generating sim model fit stats using bootstrap
  {
    all_juries_data$sim_verdict <- NA
    for(i in 1:nrow(observed.deliberations))
    {
      # simulate verdicts given observed starting conditions
      all_juries_data$sim_verdict[i] <- deliberate(g_votes=round(all_juries_data$jury_size[i]*all_juries_data$prop_jurors_g[i], 0), 
                                                          jury_n=all_juries_data$jury_size[i])
    }
    # initial model
    formal_glm = glm((sim_verdict=="G") ~ prop_jurors_g, family="binomial", data=subset(all_juries_data, !is.na(guilty_verdict)))
    formal_estimates[, ,j] <- summary(formal_glm)$coefficients
    formal_fit_stats[j, ] <- fit_stats_glm(formal_glm)
    # controlled analysis
    z_formal_glm = glm((sim_verdict=="G") ~ prop_jurors_g + death_penalty + six_person_jury, family="binomial", 
                         data=subset(all_juries_data, !is.na(guilty_verdict)))
    z_formal_estimates[, ,j] <- summary(z_formal_glm)$coefficients
    z_formal_fit_stats[j, ] <- fit_stats_glm(z_formal_glm)
  }
  initial_model <- apply(formal_estimates, c(1, 2), mean, na.rm = TRUE)
  row.names(initial_model) <- c("Constant", "Initial juror preferences"); 
  colnames(initial_model) <- c("Estimate Std.", "Error", "z value", "Pr(>|z|)")
  initial_fit_stats <- colMeans(formal_fit_stats, na.rm = T)
  names(initial_fit_stats) <- c("LR Chi2", "Pseudo R2", "CC", "PRE")
  controlled_model <- apply(z_formal_estimates, c(1, 2), mean, na.rm = TRUE)
  row.names(controlled_model) <- c("Constant", "Initial juror preferences", "Death penalty case", "Six-person jury")
  colnames(controlled_model) <- c("Estimate Std.", "Error", "z value", "Pr(>|z|)")
  controlled_fit_stats <- round(colMeans(z_formal_fit_stats, na.rm = T), 3)
  names(controlled_fit_stats) <- c("LR Chi2", "Pseudo R2", "CC", "PRE")
  
  # output summaries of results to insert into table.
  invisible(RCPA3:::headingbox("Initial Model, Empirical"))
  round(summary(initial_empirical_model)$coefficients, 2)
  round(fit_stats_glm(initial_empirical_model), 3)

  invisible(RCPA3:::headingbox("Initial Model, Formal"))
  round(initial_model, 2)
  round(initial_fit_stats, 3)
  
  invisible(RCPA3:::headingbox("Controlled Analsis, Empirical"))
  round(summary(z_empirical_model)$coefficients, 2)
  round(fit_stats_glm(z_empirical_model), 3)
  
  invisible(RCPA3:::headingbox("Controlled Analsis, Formal"))
  round(controlled_model, 2)
  round(controlled_fit_stats, 3)


  
# ---- Figure 4: Comparing Empirical Analysis and Markov Model of Jury Deliberation ----
figure4_info <- list(name="compare empirical and markov results.png", width=5, height=3)   
if(saveFiguresDesktop) { prettyPNG(figure4_info) } else { message(messageViewerFigures) }
  
  empirical_predicted = predict.glm(initial_empirical_model, newdata=data.frame(prop_jurors_g=0:12/12), 
                                    type="response", se.fit=T)
  
  par(mar=c(2.75, 3.5, 0, 0), family="serif", mfrow=c(1,1))
  yaxis_labels <- c("0.0", "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1.0")
  plot_values <- matrix(nrow = 2, byrow = T, ncol = 13, 
                        data = c(empirical_predicted$fit, get_pG_by_k(12)))
  barplot(height = plot_values, beside = T, cex.names=.8, ann = T, col = c("gray90", "gray75"), 
          legend.text=c("Empirical analysis", "Markov model"),
          args.legend=list(x=11.5, y=1, cex=.9, box.col = "white", box.lty = 0), xlab="",
          ylab="", axes=F, ylim=c(0,1.05), space = c(0, .5))
  axis(side=2, at=seq(0, 1, by=.10), labels = yaxis_labels, las=2, cex.axis=.8)
  mtext(side=1, at=1.5 + 0:12*2.5, text = 0:12, cex=.8, line = .25)
  mtext(side = 1, text = "Jurors Initially in Favor of Guilty Verdict (k)", cex = .9, line = 1.5)
  mtext(side = 2, text = "Probability of Guilty Verdict, P(G)", cex = .9, line = 2.25)
  
if(saveFiguresDesktop) { dev.off() }
  


# ---- Table 4. Illustrative Calculation of P(G) for a Specific Value of P(g) ----
p_G_by_k <- get_pG_by_k(jury_n=12)
p_k_76 <- dbinom(x = 0:jury_n, size = jury_n, prob = .76)
p_k_64 <- dbinom(x = 0:jury_n, size = jury_n, prob = .64)
table4results <- data.frame(init_k=0:12,  p_G_by_k, 
                            p_k_76, joint_probs_76 = p_G_by_k * p_k_76, 
                            p_k_64, joint_probs_64 = p_G_by_k * p_k_64)
print_table_neatly(table4results)
round(sum(p_G_by_k * p_k_76), 3)
round(sum(p_G_by_k * p_k_64), 3)


# ---- Figure 5. Increase in the Probability of a Guilty Verdict Based on Jury Pool Preferences ----
figure5_info <- list(name="increase in probability illustration replicate.png", width=3.5, height=3.5)   
if(saveFiguresDesktop) { prettyPNG(figure5_info) } else { message(messageViewerFigures) }
  
  graphics::par(mar=c(3,3,0,0), mfrow=c(1,1), omi=base::c(0.1,0.1,0.1,0.1), family="serif")
  basic.plot.grid(main="")
  pg_values_scale <- base::seq(0,1,by=.01)
  pG_values_scale <- sate::as.jury.point(pg_values_scale, jury_n=12, pstrikes=0, dstrikes=0, accuracy=.15)
  lines(x=pg_values_scale, y=pG_values_scale, col="black", lwd=1, lty=1)
  hypo_pg = .64
  actual_pg = .76
  # this shows P(g) below x-axis
  pg_placement = 0
  text(x=hypo_pg, y=pg_placement, labels = ".64", cex=.8)
  text(x=actual_pg, y=pg_placement, labels=".76", cex=.8)
  trace_arrows(sample_pg = actual_pg, col = "gray20", lty=3)
  trace_arrows(sample_pg = hypo_pg, col = "gray20", lty=3)
  # this shows corresponding P(G) values to left of y-axis
  par(xpd=T)
  horizontal_placement <- .03
  text(y=as.jury.point(actual_pg), x=horizontal_placement, label=".828", cex=.8) #
  text(y=as.jury.point(hypo_pg), x=horizontal_placement, label=".617", cex=.8) #
  par(xpd=F)
  arrows(x0 = .2, x1 = .2, y0 = as.jury.point(actual_pg), y1 = as.jury.point(hypo_pg), 
         col = "black", lty = 1, angle = 90, length = .075, code = 3)
  text(y=mean(c(as.jury.point(actual_pg), as.jury.point(hypo_pg))), x=.21,
       label=".212 difference", cex=.8, pos=4) #
  points(y=as.jury.point(actual_pg), x=actual_pg, pch=16, col="black", cex=1) #
  points(y=as.jury.point(hypo_pg), x=hypo_pg, pch=21, col="black", bg="white", cex=1) #
  legend(x=0, y=1.05, legend=base::c("Actual trial", "Hypothetical trial"),
                   col=base::c("black", "black"), pch=c(16, 1), cex=.8, pt.cex=1,
                   box.col="#FFFFFF60", bg="#FFFFFF60")
  
if(saveFiguresDesktop) { dev.off() }
  



# ---- Figure 6: Relative Uncertainty about the Value of P(G) ----
figure6_info <- list(name="comparing precision of empirical and formal models.png", width=6, height=3)   
if(saveFiguresDesktop) { prettyPNG(figure6_info) } else { message(messageViewerFigures) }
  
  # this figure takes a while to complete because some values calculates using loops
  source("compareRelativeUncertaintyFigure.R")
    
if(saveFiguresDesktop) { dev.off() }
    



# ---- Figure 7: Probabilities of Different Outcomes with Margins of Error ----
figure7_info <- list(name="prob different outcome with MOE.png", width=4.5, height=3)   
if(saveFiguresDesktop) { prettyPNG(figure7_info) } else { message(messageViewerFigures) }
      
  graph.effect.defendant(pg_actual = .76, n_actual = 240, 
                         pg_hypo = .64, n_hypo = 312, jury_n = 12)

if(saveFiguresDesktop) { dev.off() }


# ---- Figure 8: Effect of Ineffective Defense Counsel on Probability of Death Sentence ----

figure8_info <- list(name="compact display harm Porter Washington.png", width=5.5, height=1.5)   
if(saveFiguresDesktop) { prettyPNG(figure8_info) } else { message(messageViewerFigures) }

  basicCompactResultsPlot(xlab="Increase in Probability of Death Sentence")
  # showing effect CI, Porter case
  v_placement <- .080
  porter_jury_difference <- compare.jury.stats(pg_actual=.773, n_actual=761, 
                                        pg_hypo=.605, n_hypo=761, 
                                        jury_n=12, pstrikes = 10, dstrikes = 10)$difference
  segments(x0= porter_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= porter_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=porter_jury_difference[1], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=porter_jury_difference[1], y=v_placement, pos=1, label="Fla. v. Porter", font=3, cex=.9) # 
  # showing effect CI, Washington case
  washington_jury_difference <- compare.jury.stats(pg_actual=.879, n_actual=773, 
                                            pg_hypo=.911, n_hypo=773,
                                            jury_n=12, pstrikes = 10, dstrikes = 10)$difference
  v_placement <- .035
  segments(x0= washington_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= washington_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=washington_jury_difference[1], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=washington_jury_difference[1], y=v_placement, pos=1, label="Fla. v. Washington", font=3, cex=.9) # 

if(saveFiguresDesktop) { dev.off() }
  


# ---- Figure 9. Effect of Coerced Confessions on Probability of Conviction ----
figure9_info <- list(name="compact display harm Fulminante Hopkins.png", width=5.5, height=1.5)   
if(saveFiguresDesktop) { prettyPNG(figure9_info) } else { message(messageViewerFigures) }

  basicCompactResultsPlot()
  # showing effect CI, Fulminante case
  v_placement <- .080
  fulminante_jury_difference <- compare.jury.stats(pg_actual=.463, n_actual=764, 
                                           pg_hypo=.414, n_hypo=764, 
                                           jury_n=12, pstrikes = 10, dstrikes = 10)$difference
  segments(x0= fulminante_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= fulminante_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=fulminante_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=fulminante_jury_difference["Difference in P(G)"], y=v_placement, pos=1, 
       label="Ariz. v. Fulminante", font=3, cex=.9) # 
  
  # showing effect CI, Hopkins case
  hopkins_jury_difference <- compare.jury.stats(pg_actual=.858, n_actual=917, 
                                        pg_hypo=.883, n_hypo=917,
                                        jury_n=12, pstrikes = 15, dstrikes = 15)$difference
  v_placement <- .035
  segments(x0= hopkins_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= hopkins_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=hopkins_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=hopkins_jury_difference["Difference in P(G)"], y=v_placement, pos=1, 
       label="Texas v. Hopkins", font=3, cex=.9) # 

if(saveFiguresDesktop) { dev.off() }
  

# ---- Figure 10. Effect of Improper Prosecutorial Statements on Probability of Conviction ----
figure10_info <- list(name="compact display harm Chapman Lee.png", width=5.5, height=1.5)   
if(saveFiguresDesktop) { prettyPNG(figure10_info) } else { message(messageViewerFigures) }
  
  basicCompactResultsPlot()
  # showing effect CI, Chapman case
  chapman_jury_difference <- compare.jury.stats(pg_actual=.60, n_actual=99, pg_hypo=.53, n_hypo=99, jury_n=12, 
                     pstrikes = 10, dstrikes = 10)$difference
  v_placement <- .080
  segments(x0=chapman_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1=chapman_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=chapman_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=chapman_jury_difference["Difference in P(G)"], y=v_placement, pos=1, label="Cal. v. Chapman", font=3, cex=.9) # 
  
  # showing effect CI, Lee case
  lee_jury_difference <- compare.jury.stats(pg_actual=.68, n_actual=69.5, pg_hypo=.66, n_hypo=69.5, jury_n=12, 
                     pstrikes = 20, dstrikes = 20)$difference
  v_placement <- .035
  segments(x0=lee_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1=lee_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=lee_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=lee_jury_difference["Difference in P(G)"], y=v_placement, pos=1, label="Penn. v. Lee", font=3, cex=.9) # 

if(saveFiguresDesktop) { dev.off() }
  
  
# ---- Figure 11. Effect of Varying Jury Instructions on Probability of Conviction ----
figure11_info <- list(name="compact display harm Berkowitz King.png", width=5.5, height=1.5)   
if(saveFiguresDesktop) { prettyPNG(figure11_info) } else { message(messageViewerFigures) }

  basicCompactResultsPlot()
  # showing effect CI, Berkowitz case
  v_placement <- .080
  berkowitz_jury_difference <- compare.jury.stats(pg_actual=195/300, n_actual=300, pg_hypo=159/300, n_hypo=300, jury_n=12, 
                                          pstrikes = 7, dstrikes = 7)$difference
  segments(x0= berkowitz_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= berkowitz_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=berkowitz_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=berkowitz_jury_difference["Difference in P(G)"], y=v_placement, pos=1, label="Penn. v. Berkowitz", font=3, cex=.9) # 
  
  # showing effect CI, King casepg_actual <- 1-.24 # M'Naghten, text reports percent not guilty by insansity
  king_jury_difference <- compare.jury.stats(pg_actual=1-.24, n_actual=240, 
                                     pg_hypo=1-.36, n_hypo=312, jury_n=12,
                                     pstrikes=10, dstrikes=10)$difference
  v_placement <- .035
  segments(x0= king_jury_difference$`Lower 5%`,  y0 = v_placement,
           x1= king_jury_difference$`Upper 95%`, y1 = v_placement, lty=3, col="black")
  points(x=king_jury_difference["Difference in P(G)"], y=v_placement, pch=16, col="black", cex=1.1) # 
  text(x=king_jury_difference["Difference in P(G)"], y=v_placement, pos=1, label="U.S. v. King", font=3, cex=.9) # 

if(saveFiguresDesktop) { dev.off() }
  

